////////////////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                                        //
//  -----------     Ebrahim Foulaadvand and Somayyeh Belbasi, 03 August 2012     -------------------      //
//                                                                                                        //
// The programme "FreeFall" evaluates the impact velocity of freely fallen object released from a height  //
// h above the ground for two cases of constant and position dependent gravitation field g. The Euler     //
// algorithm is implemented.                                                                                             //                                                                  //
//                                                                                                        //                                                                                                         //
////////////////////////////////////////////////////////////////////////////////////////////////////////////


#include <iostream>
#include <fstream>
#include <conio>
#include <math.h>
#include <numeric>
#include <values>
#include <vector>
#include <list>
#include <algorithm>
#include <numeric>
#include<iomanip.h>

using namespace std;

int s,t,S=50;

double Y, tau=0.001, g=9.8, R=6.37e6, deltav, h, deltah=1000, v0=0., y1, y2, v1, v2;

// R= Earth radius (m), h=falling height, deltah= height increment,

main()

{

randomize();

ofstream Vimpact1 ("Vimpact constant g.plt"); // output file for the impact velocity vs height (constant g)
ofstream Vimpact2 ("Vimpact variable g.plt"); // output file for the impact velocity vs height (variable g)
ofstream delVimpact ("del Vimpact.plt");      // output file for the impact velocity deiffernec vs height


for (int s=1;s<=S;s++){

h=s*deltah; //cout<<"h= "<<h<<endl;

y2=h; v2=v0;

y1=h; v1=v0;


//------------------------ Euler Method------------------------

Y=h; t=0;

while (Y>0.01){

y2=y2 + tau*v2;

v2=v2-(tau*g)/((1+(y2/R))*(1+(y2/R)));

Y=y2;

t+=1;

}

//cout<<"t= "<<t<<endl;

Vimpact2<<h<<"  "<<fabs(v2)<<endl;


Y=h; t=0;

while (Y>0.01){

y1=y1 + tau*v1;

v1=v1-(tau*g); //cout<<"v1= "<<v1<<endl;getch();

Y=y1;

t+=1;

}

//cout<<"t= "<<t<<endl;

cout<<endl;

Vimpact1<<h<<"  "<<fabs(v1)<<endl;

deltav=fabs(v2-v1);

//cout<<"deltav= "<<deltav<<endl; cout<<endl;

delVimpact<<h<<" "<<deltav<<endl;

} // end s


cout<<"end."<<endl;getch();

return 0;
}